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I .  INTRODUCTION 


This  report  presents  the  results  of  research  supporting  the  develop¬ 
ment  and  use  of  synthetic  aperture  radar  (SAR)  as  an  instrument  for 
making  measurement  of  the  coastal  environment.  The  specific  objectives 
of  the  research  reported  here  were  (1)  to  make  ground  truth  measurements 
of  the  near  surface  drift  profile  under  waves  and  winds,  and  (2)  to  use 
numerical  modeling  of  shallow  water  wave  fields  in  the  interpretation  of 
SAR  images  of  coastal  wave  patterns. 

The  first  part  of  the  work,  that  concerning  the  near  surface  profile 
of  drift  currents,  was  conducted  as  part  of  the  MARSEN  field  experiment 
in  September  and  October,  1978.  The  MARSEN  experiment  involved  extensive 
remote  sensing  measurements  using  a  variety  of  surface  and  airborn 
radars. (1)  Needed  to  evaluate  the  accuracy  of  radar  current  measure¬ 
ments  were  insitu  techniques  which  would  describe  the  vertical  profile 
of  currents  within  a  few  meters  of  the  sea  surface.  The  research  re¬ 
ported  herein  involved  the  use  of  Lagrangran  drift  buoys  with  drag 
localizers  placed  at  depths  of  from  .5  M  to  8  M  below  the  sea  surface. 

The  buoys  emitted  a  HF  signal  and  were  tracked  using  shore  based  direc¬ 
tion  finding  (D-F)  receivers.  Triangulation  from  two  D-F  stations 
allowed  the  position  of  the  buoys  to  be  tracked  over  distances  of  up  to 
100  km  and  for  time  periods  of  up  to  250  hours. 

The  second  part  of  the  research  concerned  the  interpretation  of  SAR 
images  of  nearshore  wave  fields  for  which  the  assumption  of  wave  spectral 
homogeneity  cannot  be  made;  i.e.,  waves  on  beaches,  waves  in  the  surf 
zone,  and  wave  diffraction  around  manmade  structures.  The  interpretation 
of  SAR  images  of  such  conditions  can  be  aided  with  a  deterministic, 
rather  than  spectral,  model  for  the  underlying  wave  field.  A  finite 
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difference  deterministic  model  of  near-shore  shallow  water  wave  fields 
was  developed.  The  wave  fields  generated  by  the  model  are  spatially 
inhomogeneous  in  terms  of  wave  direction,  wave  height  and  wave  length 
and  describe  sea  surface  height,  slope,  velocity  and  acceleration  at  a 
resolution  equal  to  that  of  SAR  images.  Thus,  the  interpretation  of  SAR 
images  of  nearshore  wave  patterns  could  be  augmented  with  synthetic  SAR 
images . 


II.  DRIFT  UNDER  WAVES  AND  WIND 


The  subject  of  this  research  was  the  near  surface  profiles  of  drift 
currents  under  waves  and  winds.  While  few  field  observations  have  been 
reported  in  the  literature  on  this  subject,  several  laboratory  studies 
have  been  conducted  which  provide  some  physical  insight  into  the  dynamic 
processes  involved. 

Shemdin  1972  presented  the  results  of  laboratory  measurements  of 

(2) 

the  near-surface  drift  profile  under  combined  wind  and  wave  action. 

It  was  shown  that  a  logarithmic  boundary  layer  is  formed  near  the  air- 
sea  interface  which  increased  the  phase  speed  of  the  surface  gravity 
waves.  The  drift  was  found  to  be  3%  of  the  free  stream  wind  velocity. 
This  drift  current  increased  the  minimum  wave  phase  speed  by  100%.  For 
a  wind  of  7.9  m/sec,  waves  having  a  frequency  between  .7  and  1.6  Hz 
increased  their  phase  speed  by  .5%  to  about  10%. 

Systematic  measurements  of  drift  currents  below  the  air-water 

(3) 

interface  under  various  wind  conditions  were  reported  by  Wu  in  1975. 
Current  was  found  to  follow  a  logarithmic  profile  near  the  surface; 
however,  immediately  below  the  surface  the  current  varied  linearly  with 
depth.  It  was  found  that  the  surface  drift  was  primarily  wind-induced, 
with  Stokes  wave  drift  only  a  small  component.  However,  this  wind 
component  was  predicted  to  decrease  in  importance  with  increasing  fetch, 
while  the  wave-induced  Stokes  drift  would  increase.  Wave  drift  would 
dominate  over  wind  drift  for  fetches  of  approximately  1  km  and  greater. 

Csanady  1978,^  re-analyzing  the  data  presented  by  Wu  1975,  has 
extended  the  "law  of  the  wall,"  familiar  from  turbulent  boundary  layer 
theory,  to  turbulent  drift  flow  along  a  sharp  density  interface.  Two 
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factors  determining  the  boundary  layer  properties  were  identified. 

First,  the  effective  viscosity  in  the  laminar  sub-layer  is  greater  than 
the  molecular  viscosity  of  the  fluid  because  of  the  presence  of  wavelets 
at  the  density  interface.  These  waves  contribute  to  a  vertical  momentum 
flux  adding  to  the  effect  of  normal  molecular  viscosity.  For  values  of 
u*  (shear  velocity)  greater  than  about  1  cm/sec,  the  formula  relating 
the  effective  viscosity  (Ve)  and  the  molecular  viscosity  (V)  is  - 


where  g  is  the  acceleration  of  gravity.  Second,  the  logarithmic  (turbu¬ 
lent)  boundary  layer  indicates  a  rough  boundary,  whose  roughness  length 
scale  is  a  function  of  the  stability  of  the  water  column.  Increasing 
stability  leads  to  a  reduction  in  the  roughness  length  scale  of  the 
boundary.  The  effect  of  stability  on  the  boundary  layer  above  the 
air-sea  interface  has  been  considered  for  some  time  (i.e.,  (5)). 

The  surface  boundary  layer  formed  by  the  orbital  motion  associated 
with  surface  waves  has  not  been  as  intensively  investigated  as  for 
steady  drift  currents.  Several  wave  theories  (linear,  Stokes,  etc.)  can 
be  used  to  predict  the  vertical  velocity  profiles  through  the  main  body 
of  the  water  column.  These  present  theories  are  only  partially  success¬ 
ful  in  agreeing  with  observations.  However,  measurements  of  near  surface 
profiles  of  surface  wave  orbital  velocities  and  comparisons  with  theories 
are  very  limited. 

In  a  series  of  papers,  Shonting^ (?)(8)  presented  the  results  of 
field  observations  of  wave  particle  motions.  A  ducted-impeller  meter 
was  used  to  measure  the  Eulerian  current  at  depths  from  .5m  below  the 
wave  trough  level  to  near  the  bottom  in  20m  of  water.  There  was  a 


5 


general  exponential  decay  in  the  wave  particle  velocities  with  depth, 
however  differences  up  to  18%  were  noted  between  the  observed  values  and 
predictions  using  Stokes  5th  order  wave  theory,  linear  wave  theory  and 

stream  function  theories.  Measurements  were  made  at  10  ft.  and  25  ft. 

(0) 

depths  from  the  mean  water  level.  Forristall,  et.  al.,  made  measure¬ 
ments  of  storm  waves  and  currents  at  -13  ft.,  -23  ft.  and  -48  ft.  below 
the  mean  water  level.  Substantial  errors  occurred  (up  to  100%)  between 
predicted  values  of  the  orbital  velocity  measured  values. 

From  this  brief  review  of  the  literature,  it  appears  that  the 
kinematics  of  the  surface  layer  are  characterized  by  strong  coupling 
between  wind  effects,  waves,  wave  induced  drift  currents  and  wind  in¬ 
duced  drift  currents.  The  implication  of  these  studies  is  that  condi¬ 
tions  in  the  ocean  will  be  significantly  different  from  those  studied  in 
the  laboratory.  Major  differences  include  the  lack  of  lateral  constraints 
in  the  ocean,  the  possible  presence  of  strong  vertical  density  gradients 
and  the  existence  of  higher  and  longer  waves.  Field  measurements,  to  be 
useful  in  illucidating  the  basic  processes  taking  place,  would  have  to 
be  conducted  in  light  of  the  laboratory  results  and  the  predicted  time 
and  space  scales  of  the  processes,  with  special  attention  given  to  the 
affects  that  are  unique  to  the  field. 


1.  Field  Procedures 


The  field  site  was  the  area  of  the  North  Sea  referred  to  as  the 
German  Bight.  The  MARSEN  plan  was  to  combine  insitu  and  remote  sensing 
instruments  in  a  study  of  coastal  oceanographic  processes.  The  drift 
experiments  were  planned  for  various  parts  of  the  MARSEN  study  area. 

The  measurement  of  the  position  of  the  drift  buoys  required  two 
stations  from  which  direction  cuts  could  be  made.  These  stations  were 
on  Helgoland  (Station  A)  and  at  the  MARSEN  headquarters  on  Sylt  at 
Westerland  (Station  B) .  These  positions  are  shown  in  Figure  1. 

The  stations  were  established  when  the  position  (longitude  and 
latitude)  and  the  bearing  of  true  North  were  determined.  The  coordi¬ 
nates  of  the  stations  were  7°52.9'E.  Long.,  54°11.3'N.  Lat.  for  Station 
A  and  8°18.8'  E  Long.,  54°  55.6'  N.  Lat.  for  Station  B.  The  known 
position  of  local  landmarks  and  a  hand  held  compass  were  used  to  estab¬ 
lish  the  true  north  bearing.  The  buoys  and  direction  finding  instru¬ 
ments  were  provided  by  Telecommunications  Enterprises,  Panama  City, 

Florida^11).  Prior  use  of  the  instruments  and  the  measurement  tech- 

02) 

niques  are  presented  inv  .  The  buoys  used  in  the  study  are  shown  in 
Figure  2. 

The  buoys  used  in  the  experiments  are  listed  in  Table  1.  The  buoys 
are  distinguished  by  their  transmitting  frequency.  Buoys  were  used  for 
transmission  tests,  training  of  observers,  double  buoy  launches  and 
multiple  buoy  launches. 

The  buoys  were  deployed  with  the  drogue  located  at  various  depths 
below  the  sea  surface.  Once  the  buoys  were  launched  at  sea,  the  direc¬ 
tion  fixes  were  taken,  nominally,  at  0900,  1100,  1300,  and  1500  each  day 
until  reception  of  the  buoy  transmission  ceased.  The  buoy  measurements 
are  listed  in  Table  2. 
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2.  Analysis  of  Data 

The  locations  of  the  drift  buoys  were  determined  using  the  known 
positions  of  the  D-F  stations  and  the  measured  true  bearings  from  Station 
A  and  Station  B.  The  variables  used  in  deriving  the  position  of  the 
buoy  are  given  in  Figure  3.  The  target  bearings  from  Station  A  and  B 
are  given  as  BRA  and  BRB.  The  calculation  of  the  target  position  is 
based  upon  applying  some  relationships  involving  spherical  trigonometry, 


which  are  reviewed  first. 

For  the  spherical  triangle  shown  in  Figure  4,  the  following  rela¬ 
tionships  hold  (13): 

Law  of  Sines 

sin  AL  _  sin  BE  _  sin  GA  /jj 

sin  a  sin  b  sin  c 

where  AL,  BE  and  GA  are  angles  and  a,  b,  c  are  arc  lengths  given  in 

degrees . 

Law  of  Cosines 


cos  a  =  cos  b  cos  c  +  sin  b  sin  c  cos 

AL 

(2) 

cos  b  =  cos  a  cos  c  +  sin  a  sin  c  cos 

BE 

(3) 

cos  c  =  cos  a  cos  c  +  sin  a  sin  b  cos 

GA 

(4) 

Half-Angle  Formulas 

tan  MMS-BES)  =  -  *$$$  tan  %c 

(5) 

tan  «ALS*BES)  =  -  *“  * 

(6) 

.  i  f  _  sin^(ALS-BES)  l* 

tan  %(a  b)  -  sin\'(ALS+BE8)  1  ^ 

(7) 

The  first  step  in  the  data  analysis  was  to  compute  the  arc  length 
separation  of  the  stations  and  the  true  bearings  of  the  baseline  at  each 
station.  The  colatitudes  of  each  station  (referring  to  Figure  3)  are 


f  =  90°  -  54.18833  =  35.81670 
e  =  90°  -  54.92667  =  35.07333 
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The  difference  in  longitude  is 

DE  =  8.31333  -  7.88167  =  .43166  (8) 

and  then  the  supplement 

DES  =  180  -  .43166  =  179.56834  (9) 


Now,  using  equation  (5),  for  the  triangle  AB-NORTH  POLE  in  Figure  5, 
tan  ^(EAS-FAS)  =  -  tan  %CES  (10) 

tan  ^(EAS+FAS)  =  -  |°rji[e+fj  tan  *CES  (11) 

then  with  (EAS-FAS)  and  (EAS+FAS)  determined, 

sin  ^(EAS-FAS)  _  ,  ,  ,s  j-io-i 

tan  '  sin  % (EAS+FAS)  tan  *(e"£)  °2) 

so  that 


C  =  .77968° 


Using  an  earth  radius  of  6369.0  km,  the  separation  between  stations  was 
C  =  86.67  km 

the  true  great  circle  bearings  of  the  stations;  GCA  and  GCB  are 
GCA  =  18.5533 
GCB  =  198.5533 


which  is  found  from  the  formulas 
GCA  =  tan  *(xx)  -  tan  *(yy) 

GCB  =  tan  *(xx)  +  tan  *(yy) 

where 

_  cos  ^(LATB  -  LATA  ) 

XX  "  tan  *s(L0NG  B  -  LONG  A)  sin  ^(LATA  +  LATB) 

sin  ^(LATB  -  LATA)  _ _ 

^  '  tan  ^(LONG  B  -  LONG  A)  cos  ^(LATA  +  LATB) 

The  target  position  is  determined  using  the  known  true  bearings  of 
the  target  from  Station  A  and  B,  namely  angles  BRA  and  BRB  (See  Figure  3). 
The  steps  are  to  first  find  the  interior  angles 
AL  =  360°  +  GCA  -  BRA 


BE  =  BRB  *  GCB 
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Then,  the  values  are  found 

K1  =  sin  ^(AL  +  BE)  cos  C 
K2  =  sin  %(AL  -  BE)  cos  C 
and  zz  =  tan  *(K2/K1) 

Then,  the  angle  GA  is 

GA  =  2  cos  (  /cos  y) 

and  the  sides  a  and  b  are 

.  -1  ,sin  AL  sin 

8  *  8,11  ( - sin~GA 

,  4  -1  (sin  BE  Si.  C} 

sin  GA 

Then, 

cos  a"  =  cos  e  cos  a  +  sin  e  sin  a  cos (2  -  BRB) 

and  using  this  in 

vii  _  sin  (2  -  BRB)sin  a 

510  D  •  it 

sin  a 

Finally,  the  longitude  of  the  largest  is 
LONG  T  =  e  -  b" 
and  the  latitude  is 
LAT  T  =  90  -  a" 

The  final  step  was  to  compute  the  separation  between  the  position  of  the 
target  at  two  different  times  T1  and  T2,  the  bearing  of  the  movement  and 
the  average  speed.  The  separation  and  movement  bearings  are  calculated 
using  the  functions  used  to  determine  station  separation  and  bearing  of 
the  baseline.  The  average  speed  was  calculated  from 
VA  =  CX/T2-T1 

where  CX  is  the  distance  the  buoy  moved  in  the  time  interval  T2-T1. 
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3.  Results 

The  results  of  the  data  analyses  for  selected  buoy  experiments  are 
presented  in  Tables  3-6  and  Figures  6  and  7.  The  tables  give  the  position 
of  a  buoy,  the  angle  GA  of  each  measurement,  and  for  pairs  of  measurements 
the  distance  between  positions  (DIS),  direction  of  movement  (DIR)  and 
the  average  speed  (SPD).  The  position  of  Buoy  396  for  the  period  up  to 
112  hours  after  launch  is  shown  in  Figure  6. 

The  buoy  experiment  shown  in  Figures  6  and  7  involves  2  buoys  (#396 
and  #397)  for  which  the  drogue  depth  was  .5  m  and  8  m,  respectively. 

The  Tables  4  and  5  indicate  that  during  low  wind  speed  when  currents 
were  dominated  by  tides  (hours  from  35  to  86)  the  ratio  of  the  shallow 
drogue  speed  to  the  deep  drogue  speed  varied  both  above  and  below  the 
value  1,  with  an  average  value  of  .86.  However,  during  the  strong  wind 
period  (hours  from  133  to  169)  the  ratio  was  consistently  above  1.0  and 
averaged  1.75.  This  trend  can  be  seen  in  Figure  7,  where  the  shallow 
buoy  396  position  changes  are  consistently  larger  than  the  deeper  buoy 
397. 

The  accuracy  of  the  profile  measurements  depend  upon  the  accuracy 
of  the  measurement  of  buoy  position.  For  a  bearing  accuracy  of  1°  , 
the  position  accuracy  at  50  km  is  .87  km  .  Most  of  the  distances 
between  buoy  fixes  were  of  the  order  10  km  ,  although  some  separations 
were  as  great  as  39  km.  Thus*  the  accuracy  of  the  speed  measurements 
were  of  the  order  of  10  7.  .  A  10  7.  error  in  the  speed  measurements 
translates  to  about  a  20  7.  error  in  the  ratio  of  the  speeds.  Therefore 
the  ratio  measured  for  the  buoy  pair  cited  above  is  .86  -  .17  for 
light  winds  and  1.75  -  .35  for  strong  winds. 


11 


The  speed  profile  of  drift  under  waves  and  winds  has  been  described 
for  a  PM  spectrum  as  decreasing  as  exp(2kz),  where  k  is  the  wave  number 
and  z  is  the  depth  (negative)  below  the  surface  (18).  For  the  ratio 
of  speeds  between  two  depths  z^  and  z^  ,  the  formula  becomes  exp^z^-z^) )• 
For  z^-,5  m  and  z^-S  m  and  a  wave  frequency  of  .16  Hz,  the  predicted 
ratio  is  2.16  .  The  observed  ratio  under  high  winds  was  1.75  -  .35  . 
Therefore  the  drift  measurements  tend  to  support  the  exponential 
decay  of  speed  below  the  surface.  Tidal  currents  ,  which  were  ignored 
in  the  theoretical  formula  ,  would  tend  to  reduce  the  actual  ratio 
below  the  theoretical  value. 


III.  SYNTHETIC  SAR  IMAGES 


Synthetic  aperture  radar  (  SAK  )  is  now  an  operational  system  on 
Marine  Corps  reconnaisance  aircraft  in  the  form  of  the  Goodyear  AN/UPD-4  (14). 
Many  environmental  parameters  of  importance  to  Naval  operations  can  be 
sensed  by  this  system.  In  particular  the  system  is  able  to  image  nearshore 
waves  and  surf,  parameters  which  are  needed  to  define  the  environmental 
conditions  of  many  Naval  Inshore  Warfare  operations.  While  considerable 
effort  is  being  devote  to  develop  the  capability  of  SAK  to  measure  deep¬ 
water  waves,  relatively  little  research  is  being  directed  toward  the  special 
problems  of  measuring  nearshore  waves  and  surf.  The  purpose  of  this  study 
was  to  develop  a  method  for  deterministically  analysing  SAK  images  of 
nearshore  waves  and  surf. 

The  use  of  SAR  to  determine  wave  properties  in  deep  water  is  based 
upon  describing  the  sea  surface  with  a  directional  spectrum  of  waves  (15,16). 
Typically  an  area  of  the  image  is  Fourier  transformed  and  compared  with  the 
ground  truth  directional  wave  spectrum.  This  approach  requires  that  the  waves 
over  the  area  of  the  portion  of  the  image  analysed  have  a  homogeneous 
two-dimensional  spectrum.  However  wave  patterns  in  shallow  water  are  not 
homogeneous  due  to  the  effects  of  such  processes  as  refraction,  shoaling 
and  attenuation.  The  approach  taken  in  this  study  was  to  develop  ground 
truth  for  an  area  in  shallow  water  by  constructing  the  sea  surface  topog¬ 
raphy  and  kinematics  using  a  computer  model.  This  information  could  then 
be  used  ,  using  various  radar  models,  to  construct  synthetic  SAR  images 
which  could  be  compared  with  the  actual  images.  This  approach  has  the 
advantage  that  it  can  be  used  around  man-made  objects  and  structures 
which  are  so  prevalent  in  the  coastal  areas  of  the  world  • 
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1.  Wave  Model  Approach 

All  of  the  classical  analytical  wave  theories  have  the  common  limi¬ 
tation  of  referring  to  simple  geometric  situations,  i.e.,  a  plane  long- 
crested  wave  moving  over  a  flat  or  constant  sloping  bottom.  The  general 
non-linear  equations  of  motion  with  an  arbitrary  bottom  shape  are  not 
soluble  analytically.  In  order  to  treat  the  problem  of  real  waves  in 
the  detail  required  to  compare  with  SAR  images,  we  are  compelled  to  intro¬ 
duce  discretization  of  the  space  ans  time  dimensions  ,  and  solve  for  the 
dependent  variables  at  all  the  necessary  points  of  the  four-dimensional 
space-time  grid.  This  means  converting  the  analytic  equations  to  finite 
differnce  form  and  solving  these  numerically. 

However,  the  numerical  problem  posed  by  the  general  finite  difference 
equations  is  still  formidable  on  account  of  its  spatially  three  dimensional 
character.  In  the  case  of  a  homogeneous  fluid  which  is  isothermal  and  in¬ 
compressible  there  are  four  dependent  variables:  u,v,w,p  --  representing 
the  three  components  of  the  velocity  and  the  pressure,  and  the  water 
level  .  These  must  be  calculated  for  every  time  step  which  over  a 
grid  of  say  100  by  100  points  ,  will  require  excessive  computer  storage 
and  computational  time. 

A  possible  way  out  of  these  difficulties  is  to  use  the  general  charac¬ 
teristics  of  the  vertical  variation  in  the  dependent  variables  and  inte¬ 
grate  the  equations  of  motion  from  the  bottom  to  the  surface.  This  removes 
one  of  the  spatial  variables  from  the  problem. 
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2.  Derivation  of  the  Vertically  Integrated  Equation*  for  the  General 
Wave  Problem 

The  three-dimensional,  ensemble-averaged  equations  of  motion  for  a 
homogeneous,  incompressible  fluid,  with  turbulent  stresses  on  vertical 
planes  and  molecular  viscous  stresses  neglected,  can  be  written 


3u  +  3v  +  3w  *  0 
3x  3y  3z 


3u 

+  3u^ 

+ 

3uv 

+ 

3uw  - 

3t 

3x 

3y 

3z 

3v 

3vu 

+ 

3v2 

+ 

3vw  + 

3t 

3x 

3y 

3z 

3w 

.  3wu 

i 

3wv 

A. 

3w2 

3t 

+ 

*T 

3y 

T 

df  + 

-1  3p  -  1  3t«z 
p  3x  p  3z 


P  3y  p  3z 


1  3^zz 

R  -  T  TT- - 


(1) 


where  is  the  earth's  angular  velocity  vector,  x  ,  x  ,  x  are 

xz  yz  zz 

turbulent  stress  components  on  the  horizontal  plane  in  the  directions  x, 


y,  z,  respectively,  g  is  the  acceleration  due  to  gravity;  and  other 
symbols  have  been  defined.  The  Coriolis  components  of  acceleration  have 
been  simplified  by  omission  of  terms  involving  w,  on  the  assumption  that 


w  <<  u,  v. 

We  eliminate  the  third  momentum  equation  at  the  outset  by  adopting 
the  pressure  distribution  from  small  amplitude,  frictionless  wave  theory 

(see,  for  example  Ippen  (1966)): 

p  *  -gpz  +  Rpcosh[k(h  +  z)].c  (2) 

cosh[kh] 


In  the  derivation  of  equation  (2)  the  surface  boundary  condition 
was  applied  not  at  the  actual  surface  z  ■  C,  but  at  z  *  0.  Conse¬ 
quently,  the  equation  has  validity  only  from  z  *  -h  to  z  ■  0.  For 

positive  z  one  must  assume  a  hydrostatic  pressure  distribution. 

Figure  8  illustrates  the  definitions  of  z,  h  and 
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Write  ^(z)  ,  coah[k(h  4  z)1  -h  £  z  £  0 

cosh[kh] 


(3) 


tp  has  been  termed  the  "pressure  response  factor".  From  equations  (2) 
and  (3 ) 


1  Ip  -  »P  *1  .  I  M  “  g#  15. 

p  3x  3x  p  3y  3y 


With  the  last  two  relations,  and  replacing  2H  by  the  Coriolis 

z 

parameter  f,  the  equations  to  be  solved  become 

3u  +  3v  +  3w 
3x  3y  3z 

3u  +  3u2  +  3uv  +  3uw  -  fv  ■  -g\J>  3g  -  1  3^xz  ^ 

3t  3x  3y  3z  3x  p  3z 

3v  +  3vu  +  3v2  +  3vw  +  fu  ■  -g^  3^-1  3*Eyz 

3t  3x  3y  3z  3y  p  3z 


Equations  (4)  must  be  integrated  over  z  between  the  limits  -  h  and 
£.  Consider  firstly  the  integration  of  Here,  in  accordance  with  the 
previous  remarks,  the  upper  limit  must  be  changed  to  z  *  0.  The  error 
will  not  be  serious  if  the  small  amplitude  condition:  |^|«  h  ,  holds. 
From  equation  (3)  0 

/  <Mz  -  1  tanh[kh]  (5) 

-h  k 


Now,  another  result  of  the  linear,  or  small  amplitude  theory  is 


that 


a2  *  gk  tanh[khj 


(6) 


where  <j  is  the  angular  frequency  of  the  wave,  a  constant.  Eliminating 
tanh  [kh]  from  equation  (5)  using  equation  (6),  we  then  obtain  approxi- 

“tely  /C  *  95  d,  -  95 

-h  ‘5x  k*  ^ 


(7) 
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By  definition  of  the  phase  speed,  S 

o2  -  S2 


For  long  waves,  when  *l>  ■  1  is  a  good  approximation,  the  right  hand 

side  of  equation  (7)  is  replaced  by 

g(C  +  h)  li 
3x 

Boundary  conditions  that  must  be  employed  in  the  vertical 


integrations  are: 


3£  +  uH  +  v9£  -  w 
3t  3x  3y 


z  *  Z 


3h  +  u3h  +  v  3h^  *  -w 
3t  3x  3y 


Using  Leibnitz's  Rule  with  conditions  (9)  and  (10)  and  relations  of 
type  (7),  we  vertically  integrate  equations  (4)  between  the  limits  z  *  -h 


and  z  «  £  to  yield 

^  +  3_  udz  +  3_  vdz  “  0 
3t  3x  -h  3y  -h 


3_  I1  udz  +  u2dz  +  3_  uvdz  -  f/’vdz  -  -cr2  3£  +  T®  -  Tx 

3t  -h  3x  -h  3y  -h  -h  k7  - - - 


3x  -h 


3y  -h 


3  vdz  +  ^ .  vudz  +  |—  ^  v2dz  +  f/**  udz  *  -a2  3£  +  T®  - 

37  -h  wX'h  9y-h  -h  i?  37  -* - 2. 
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Tx  and  tx  are  the  x-components  of  the  surface  and  bottom  stresses. 


respectively.  Similarly  with  the  y-components. 
We  now  follow  Leendertse  (17)  t  in  writing 
u  -  U(1  +  eu) 
v  •  V(1  +  Cj 


(12) 


where 


V 


(13) 


U  “  1  udz  , 

h  +  C  -h 


vdz 


h  +  c  -h 

Thus  U  (x,  y)  and  V  (x,  y)  are  the  vertically  averaged  values  of  u 
(x,  y,  z)  and  v  (x,  y,  z),  respectively;  and  (x,  y,  z)  and  ev(x,  y, 
z)  express  the  fractional  deviations  of  u  and  v  from  U  and  V,  respec¬ 
tively,  at  x,  y,  z. 

Define  the  total  height  of  the  water  column 

H  -  h  +  X. 


Then,  using  equations  (12)  and  (13),  the  left  hand  sides  of  equa¬ 
tions  (11)  become 

+  3HU  +  3HV 
3t  3x  3y 

3HU  +  3_  {U2(H  +  /?e 2  dz)}  +  g_  {uv(h  +  f*e  f  dz)}  -  fHV 
3t  dx  -h  3y  -h 

3HV  +  {VU(H  +  /i’euevdz) }+  i_  {V2(H  +  /Ce2dz)  }  +  fHU 
3t  3x  3y  -h  V 


Assume  isotropy  of  the  horizontal  velocity  stratification. 

e  (x,  y,  z)  ■  g.  (x,  y,  z)  -  E(x,  y,  z) 
u  v 


and  let 

H  +  /?  £2dz  -  aH 
-h 


That  i6 

(14) 

(15) 


a(x»  y*  t)  may  be  termed  the  "local  velocity  stratification  factor" 
at  time  t.  Always 

a  >  1 


With  assumption  (14)  and  definition  (15),  and  after  some  reduction 
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the  equations  governing  the  flow  become 


3£  +  3HU  +  3HV  -  0 
3t  3x  3y 


(16) 


3U  +  aU  3U  +  aV3U 
3t  l)x  3y 


fV 


-21  3C  +  t“ 

Ph  15x 


-  Tb  -  {u23a  +  UV  3a  +  U 
- *  H 


(1  -  a 


PH 


3;} 

It 


3V  +  aU  3V  +  aV3V  +  fu  -  -a2  3£  +  ts  -  rb  -  {VU  3a  +  V2  3a  +  V  (1  -  a)3?) 

3t  3x  3y  k^H  3y  _y _ 3x  3y  H  3t 

PH 


2  2 

If  in  equations  (16)  we  put  a  ■  1  and  replace  o  /k  H  by  g,  we 
obtain  the  Ha  nsen  long  wave  equations.  From  (6)  we  see  that  the  second 
replacement  amounts  to  the  equality 

tanh  [kh]  -  kH 

This  will  be  closely  approached  for  shallow  water  (or  long)  waves 
provided  |?1  <<  h. 

3.  Evaluation  of  the  Velocity  Stratification  Factor 

It  is  desired  to  evaluate  a  in  closed  form,  using  the  results  of 
linear  wave  theory. 

On  the  same  level  of  approximation  as  equation  (2) 

u  *  gk  cosh[k(fi  +  z)] .  C  -h  z  <_  0  (17) 

o  cosh[kh] 

From  equations  (3),  (5)  and  (17)  we  have 

/\idz  *  g£tanh[kh] 

-h  0 

Hence,  approximately,  by  definition  (13)  and  the  last  result 

U  *  g£tank[kh]  (18) 

OH 

and  therefore,  approximately 

C  »  UoH 

gtanh[kh] 

Substituting  the  last  approximation  for  C  in  (17)  yields 

u  »  UkH  cosh[k(h  +  z) 1 
sinhfkhj 
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whence,  comparing  with  the  first  of  equations  (12) 

1  +  e  «*  kH  cosh[k(h  +  z) ]  (19) 

sinh[kh] 

Now,  by  definition  of  £ ,  we  have 

(1  +  e)2  dz  «  H  +  e2  dz 
-h  -h 

Hence,  from  equations  (15),  (19),  and  the  last  result 

aH  =  k2H2  cosh2[k(h  4-  z)]  dz 
-h  sinh2[kh] 

from  which  2 

a  -  1  (kH)  (1  +  sinhC2kH^ 

2  (sinh[kh])2  2  kH 

It  will  be  desirable,  for  computational  reasons,  that  ot  should  be  a 
function  of  position  only.  To  remove  the  time  dependency  of  ot ,  we 

substitute  h  for  H  in  the  last  result,  and  so  obtain 

a  -  _1  (kh)2  ,  (1  +  sinhfekh-])  (20 

2  (sinh[kh])  2  kh 

Equation  (20)  again  only  holds  with  the  small  amplitude  assumption. 
Some  computed  values  of  a  between  the  asymptotic  limits  1  and  kh/2 
are  shown  in  Table  8 . 

Table  8 

Some  Computed  Values  of  a 


kh 

a 

Long  Waves 

0 

1 

(Shallow  water) 

0.3 

1.00 

Intermediate 

1 

1.02 

Waves 

2 

1.18 

3 

1.53 

It  appears  that  the  Hansen  equations  are  adequate  well  into  the 
"intermediate  wave"  zone,  and  should  yield  fair  accuracy  for  wavelengths 
as  short  as  3h. 

4.  The  Two-Dimensional  Numerical  Model 


Prior  to  the  present  study  a  two-dimensional,  vertically-averaged, 
long  wave  computer  model  had  been  developed  by  McHugh  (18  1,  utilizing 
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a  mixed  implicit-explicit  method  of  solving  the  finite-difference 
equations  derived  from  Leendertse  (1967).  The  discretization  was  that 
shown  in  Figure  9. 

A  full  account  of  the  genesis  of  the  model,  which  was  designed  to 
accomodate  stretching  of  the  grid  lines  and  a  time-varying  land/water 
boundary,  is  given  in  McHugh  (1976).  The  model  was  subsequently  modi¬ 
fied  to  allow  randomly  spaced  grid  lines,  and  this  is  the  form  utilized 
in  the  present  work. 

In  order  to  generalize  the  long  wave  model  to  simulate  waves  of 
arbitrary  length,  it  was  merely  necessary  to  insert  into  the  momentum 
equations  the  factors  a  and  the  correction  terms  in  braces  as  shown  in 

equations  (16);  also  to  replace  g  in  the  pressure  driving  terms  by 

2  2  2 
cT/k  H  -  S  /H. 

The  bottom  stress  was  represented  through  the  Ch&zy  coefficient, 

Cu.-j.t').  T V.*.; 

=  0KU(U2  +  V2)1/2  .  Tb  »  pgV(U2  +  V2)1/2 

x  c2  y  — & - 

C  values  were  computed  from  Manning’s  formula: 

c  ■  1-486  H1/6  (H  in  ft  ) 

n 

The  resultant  finite  difference  equations,  which  are  very  involved, 
have  been  presented  as  Appendix  1.  Referring  to  this  appendix,  the 
superscript  n  (1,  2,  3...)  denotes  the  time  level  nAt.  Barred  quan¬ 
tities  with  subscript  u  or  v  denote  two-point  interpolated  values  at  the 
positions  of  U  or  V  indicated.  Quantities  V  *  .  and  U  *  .  denote  four- 

point  Interpolations  of  V  and  U,  respectively,  centered  at  the  locations 

of  U  .  and  V  . .  All  of  these  interpolations  take  Into  account  the 
itj  t,j 

possible  non-uniform  grid  spacing. 
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The  quantities  I|,  J*  etc.  are  aultlpllers  for  determining  deriv¬ 
atives  when  the  grid  spacing  is  irregular  or  stretched.  In  equation 
A(i)  for  instance,  for  irregular  spacing: 

r,  1 


T  » 

i+1/2 


4+1  -  x. 


The  notation  derives  from  the  form  when  1  •  I(x)  and  J  ■  J(y)  are 

stretching  functions  such  that 

9__3_dI.3_.a_dJ 
3x  31  dx  3y  3J  dy 

or,  in  finite  difference  operator  form,  writing  I*  for  dl/dx  and  J'  for 


dJ/dy 


A_«=  A_  1’ 
Ax  Al 


A_  J* 
AJ 


Lower  case  subscripts  have  been  used  in  the  equations  for 
convenience . 

Quantities  S.  .  and  a  were  calculated  at  the  centers  of  the 
J 

corresponding  grid  squares.  Two  sorts  of  derivatives  of  ot  were  used: 
those  centered  on  the  ^-points,  and  those  centered  on  the  U  and  V-points. 
The  latter  are  distinguished  by  additional  subscripts  u  or  v  in  the 
manner  already  discussed.  Thus  for  example 
“xi.j  *  _ ai+l,J  ~  al-l,J 

l/2(xi+2  +  xi+1  -  xt  -  x^ 

^XUi.j  *  «i,J  ~  ai-l,j 

l/2(x1+1  -  xi_i) 

with  analogous  expressions  for  oiy  and  Oyv  . 

i»j  i*J 

For  the  constant  boundary  model,  k(x,  y) ,  S  (x,  y) ,  and  a(x,  y) 
are  calculated  once  and  for  all  time  over  the  entire  wet  field  from 
equations  (6),  (8)  and  (20),  respectively.  Newton's  iterative  method 
was  found  adequate  for  extracting  k  from  equation  (6).  The  first  guess 
of  k  was  the  deep  water  value,  k^,  given  as  input  to  the  computer  program. 
When  the  wet  boundary  changes  with  time  (as  for  Instance  on  a  beach) 
there  is  the  necessity  of  defining  S  and  a  in  the  newly  flooded  region. 

The  strategy  in  the  present  study  was  to  assume  a  default  value  o  ■  1, 

2 

and  a  default  value  S  /H  •  g  in  the  initially  dry  regions  (h  <  0). 
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The  implicit-explicit,  double  sweep  solution  technique  ha6  been 
described  by  Leendertse  (1967)  and  McHugh  (1976).  Its  principal  stages 
are:  (i)  Solve  implicitly  for  +  *  and  Un  +  *  ^long  each  grid  row. 

(ii)  Solve  explicitly  for  Vn  +  *  in  each  square  of  the  wet  field,  (iii) 

Solve  implicitly  for  £n  +  ^  and  Vn  +  ^  along  each  grid  column,  (iv) 

Solve  explicitly  for  Un  +  ^  in  each  square  of  the  wet  field.  Repeat  the 
procedure  for  the  next  pair  of  time  steps,  and  so  on.  If  the  land/water 
boundary  varies  in  time,  this  must  be  established  periodically — eg., 
every  two  time  steps — and  allowance  made  for  volume  transport  into 
previously  dry  areas  (see  Appendix  2) . 

A  principle  of  the  computation  is  to  ignore  derivatives  when  they 
cannot  be  computed  from  the  central-difference  formula  owing  to  the 
proximity  of  a  land  boundary.  This  is  found  (Leendertse  (1967))  to  be 
beneficial  for  numerical  stability. 

Finally,  the  necessary  boundary  conditions  required  on  a  row  or 
column  are  either  the  water  level  C,  or  the  input/output  velocity  U  or  V 
at  all  times.  At  a  land  boundary  U  ■  V  ■  0. 

At  the  start  of  the  model  it  is  customary  to  assume  S  ■  O.cv'J  O 

everywhere  inside  the  wet  field. 

5.  Simulation  for  Rona  Rock 

The  numerical  model  was  applied  to  a  nearshore  area  on  the  leeward 
side  of  Rona  Rock,  a  small  island  imaged  using  the  ERIK  SAK.  The  image 
was  taken  of  waves  approaching  the  island  whose  height  was  5.5  m  and  whose 
deepwater  wavelength  was  211  m  (  R.  Schuchman  ,  personal  communication  ). 
The  leeward  side  of  the  island  was  modeled  using  the  bathymetry  shown 
in  Figure  10.  The  wave  period  was  taken  to  be  11.6  sec.  The  grid 
spacing  used  in  the  model  was  25  m  which  for  a  grid  of  120  grid  spaces 
means  that  the  model  grid  covers  an  area  3  km  square.  Because  of  the 

I 
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importance  of  several  wave  kinematic  parameters  to  the  SAR  imaging  mechan¬ 
ism,  several  variables  were  computed  at  each  grid  point.  These  parameters 
were;  sea  surface  elevation,  horizontal  wave  orbital  speed  in  the  x  and  y 
directions,  the  sea  surface  slope  in  the  x  and  y  directions  and  the 
horizontal  acceleration  and  vertical  acceleration.  The  input  wave  was  assumed 
to  be  long  crested  in  deep  water. 

6.  Results 

The  results  of  the  computation  are  given  in  Appendix  3.  As  an  example 
of  the  output  from  the  model  the  sea  surface  elevations  are  shown  in 
Figure  11.  The  Figure  shows  elevations  above  mean  sea  level  as  dark 
grid  points.  The  refraction/diffraction  patttern  of  the  waves  in  the 
lee  of  the  island  can  be  seen.  If  the  SAR  images  were  solely  a  function 
of  sea  surface  topography,  then  the  SAR  image  would  be  directly  comparable 
to  the  model  output. 
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TABLE  1.  BUOY  RECORD 


Buoy 

Drag 

Buoy 

Launch^ 

2 

Launch 

Number  Depth 
(M) 

Frequency 

(KHZ) 

Time/Date 

Position 

389 

4 

4146.5 

0615/8  Oct. 

7.3167E,  54.9167N 

391 

4 

4130.0 

0615/8  Oct. 

7.3167E,  54.9167N 

392 

2 

4229.4 

Test  Buoy 

393 

2 

4152.9 

Test  Buoy 

394 

3 

4193.3 

Test  Buoy 

395 

8 

4234.0 

0615/8  Oct. 

7.3167E,  54.9167N 

396 

.5 

4183.7 

1905/26  Sept. 

7.5035E,  54.4975N 

397 

8 

4217.7 

1905/26  Sept. 

7.5035E,  54.4975N 

398 

4 

4222.8 

0900/22  Sept. 

7.7500E,  54.7500W 

399 

1 

4117.0 

0900/22  Sept. 

7.7500E,  54.7500W 

400 

1 

4100.4 

0900/22  Sept. 

7.7500E,  54.7500N 

401 

2 

4288.5 

0900/4  Oct. 

7.1833E,  54.6833N 

402 

2 

4187.6 

0615/8  Oct. 

7.3167E,  54.9167N 

403 

2 

4286.6 

0615/8  Oct. 

7.3167E,  54.9167N 

408 

5.3 

4174.1 

Test  Buoy 

409 

3 

4280.2 

0900/22  Sept. 

7.7500E,  54.7500N 

Time  is  local  experiment  time,  which  was  European  Standard  Time. 

Launch  position  given  in  degrees  of  east  longitude  and  north 
latitude.  Values  given  in  degrees  and  decimal;  i.e.,  7°30'  =  7.5000° 
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TABLE  3.  RESULTS  OF  OBSERVATIONS  FOR  BUOY  394 
(Drogue  Depth  3M) 


Hours  from 

Position 

GA 

DIS 

DIR 

SPD 

Launch 

Long.,  Lat. 

(KM) 

(CM/S) 

14.25 

7.9559,54.787 

128 

4.46 

222 

145.7 

15.1 

7.9093,54.7576 

127 

1.11 

279 

30.92 

16.1 

7.8920,54.7593 

125 

5.61 

350 

6.74 

39.2 

7.8772,54.8090 

115 

3.50 

174 

4.88 

59.1 

7.8225,54.8089 

109 

2.93 

174 

35.4 

61.4 

7.8271,54.7826 

114 

2.03 

64 

33.3 

63.1 

7.8556,54.7906 

115 

7.78 

359 

10.8 

83.2 

7.8537,54.8606 

102 

7.24 

255 

105.8 

85.1 

7.7442,54.8440 

97 

7.47 

55 

98.8 

87.2 

7.8398,54.8824 

97 

30 


TABLE  4.  RESULTS  OF  OBSERVATIONS  FOR  BUOY  396 
(Drogue  Depth  ,5M) 


Hours  from 

Position 

GA 

DIS 

DIR 

SPD 

Launch 

Long.,  Lat. 

(KM) 

(CM/S) 

11.6 

7.5810,54.4440 

104 

9.45 

151 

107.2 

14.2 

7.6519,54.3702 

108 

4.34 

345 

6.0 

34.2 

7.6348,54.4080 

109 

4.16 

324 

55.0 

36.3 

7.5973,54.4385 

106 

5.96 

137 

85.0 

38.2 

7.6597,54.3990 

112 

3.78 

124 

48.9 

40.4 

7.7083,54.3802 

119 

10.27 

319 

15.9 

58.4 

7.6042,54.4501 

107 

2.35 

236 

32.6 

60.35 

7.5740,54.4382 

103 

2.40 

23 

36.1 

62.2 

7.5886,54.4580 

105 

3.71 

164 

554.3 

64.1 

7.6046,54.426 

106 

4.67 

305 

6.3 

84.2 

7.5453,54.4501 

100 

7.34 

122 

47.0 

88.45 

7.6408,54.4166 

110 

7.10 

2.5 

11.2 

106.1 

7.6457,54.4780 

113 

1.85 

333 

22.4 

108.3 

7.6332,54.4934 

112 

8.40 

330 

55.6 

112.5 

7.5680,54.5590 

104 

14.9 

332 

21.0 

132.3 

7.4592,54.6772 

90 

3.76 

341 

52.2 

134.2 

7.4406,54.7093 

87 

4.39 

167 

61.1 

136.2 

7.4557,54.6707 

90 

15.06 

352 

21.1 

156 

7.4255,54.8051 

80 

11.06 

285 

76.8 

160 

7.2587,54.8311 

69.5 

9.67 

39 

14.9 

178 

7.3541,54.8986 

69.5 
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TABLE  5.  RESULTS  OF  OBSERVATIONS  FOR  BUOY  39? 
(Drogue  Depth  8M) 


Hours  from 

Position 

GA 

DIS 

DIR 

SPD 

Launch 

Long.,  Lat. 

(KM) 

(CM/S) 

10.2 

7.5911,54.4497 

105 

6.86 

157 

48.9 

14.1 

7.6328,54.3930 

108 

1.39 

83 

20.4 

16.0 

7.6541,54.3945 

111 

7.51 

327 

11.6 

34.0 

7.5907,54.4511 

105 

.51 

192 

7.1 

36.0 

7.5890,54.4466 

105 

4.09 

131 

56.8 

38.0 

7.6364,54.4223 

110 

5.14 

147 

71.3 

40.0 

7.6792,54.3834 

114 

8.25 

320 

12.7 

58.0 

7.5980,54.4407 

106 

4.23 

211 

58.7 

60 

7.5641,54.4081 

100 

6.28 

13 

87.3 

62 

7.5864,54.4631 

105 

3.62 

142 

50.3 

64 

7.6203,54.4372 

109 

8.19 

295 

11.4 

84 

7.5056,54.4689 

96 

8.17 

121 

56.7 

88 

7.6132,54.4304 

107 

10.49 

11 

14.6 

108 

7.6450,54.5229 

113 

6.11 

305 

84.9 

110 

7.5676,54.5548 

104 

20.17 

339 

25.5 

132 

7.4573,54.7247 

87 

3.37 

274 

46.9 

134 

7.4047,54.7269 

84 

2.28 

25 

31.7 

136 

7.4195,54.7455 

84 

9.47 

321 

13.1 

156 

7.3270,54.8120 

74 

4.47 

42 

31.1 

160 

7.3736,54.8419 

74 

5.98 

328 

9.2 

178 

7.3240,54.8876 

69 

4.04 

220 

56.2 

180 

7.2832,54.8598 

69 

32 


TABLE  5  Continued 


Hours  from 
Launch 

Position 

Long.,  Lat. 

GA 

DIS 

(KM) 

DIR 

SPD 

(CM/S) 

202 

7.1176,55.0048 

55 

19.3 

327 

24.3 

250 

6.8530,55.3230 

37 

39.15 

334 

22.6 

253 

6.6714,55.4111 

33 

15.07 

311 

139.56 

TABLE  6.  RESULTS  OF  OBSERVATIONS  FOR  BUOY  408 
(Drogue  Depth  5.3M) 


Hours  from 

Position 

GA 

DIS 

DIR 

SPD 

Launch 

Long. ,  Lat. 

(KM) 

(CM/S) 

14.2 

7.8609,54.7259 

126 

3.06 

43 

85.1 

15.2 

7.8937,54.7458 

127 

4.07 

12 

118.6 

16.3 

7.9085,54.7872 

122 

4.52 

6 

7.01 

34.2 

7.9157,54.8276 

115 

5.03 

226 

34.1 

38.3 

7.8590,54.7962 

115 

3.73 

314 

4.9 

59.3 

7.8174,54.8198 

107 

2.09 

163 

27.6 

61.4 

7.8270,54.8018 

111 

2.62 

150 

38.3 

63.3 

7.8471,54.7814 

116 

9.94 

17 

13.8 

83.2 

7.8946,54.8666 

104 

7.57 

251 

105.2 

85.2 

7.7824,54.8445 

100 

10.09 

63 

127.4 

87.4 

7.9224,54.8862 

102 
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TABLE  7.  COMPARISON  OF  DRIFT  SPEEDS 


s  from 

Speeds 

(CM/S) 

Ratio  of  S 

ch 

396 

397 

396/397 

35 

55 

7.1 

7.7 

37 

85 

56.8 

1.5 

39 

48.9 

71.3 

.68 

49 

15.9 

12.7 

1.25 

59 

32.6 

58.7 

.56 

61 

36.1 

87.3 

.41 

63 

54.3 

50.3 

1.08 

74 

6.3 

11.4 

.55 

86 

47.0 

56.7 

.83 

133 

52.2 

46.9 

1.11 

135 

61.1 

31.7 

1.93 

146 

21.1 

13.1 

1.61 

158 

76.8 

31.1 

2.47 

169 

14.9 

9.2 

1.62 

VI.  FIGURES 


1, 


Figure  1.  Field  site  and  D-F  stations. 


DROGUE  DEPTH 


Figure  2.  Buoy  types  used  in  the  experiments 
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GRID  NUMBER 


120 


Figure  11.  Elevations  predicted  for  leeward  side  of  Rona  Rock  using 

the  numerical  model.  Dark  is  positive  and  white  is  negative 
elevation. 
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Appendix  1  -  The  Finite-Difference  Equations 
of  the  General  Vave  Model 

(See  Section  4  for  a  description  of  symbols) 


C  ■  +  <Hui+l,j  Uitl,4  •  Hu"  -  U”+^  :^1/2 


,5  i.J  i,J 


At 


+  (Hvi,j+l  Vi,j+1  ~  **vi,j  Vi,j)  Jj+l/2 


A(i) 


Un+1  -  [Un  1 

ui,j  lui,jJ  -  n+1  n  t|n 

(XU.,  .  U .  .  (Uj  . ,  ,  -  u 


At 


|)  1  .  «...  It**1  Jy® 


+  "‘.J  i  +  Vj  <"i,j« '  "i.j- 


f  v*°j  +  sj.j. .  -  Ci1.)'  +  <j 

2h<.3 


Ci.i'i  U 


+  f  +  <Cj»  'fC/  +  <v‘i?j>2>1/2  (5ri.j>2  -  K.) 


1}  Jj‘+l/2 


A(ii) 
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Vtv+^  _  fyn  1  +  av  u*n+  (V11  -  Vn  1  I  ’  ,  +  av  (vn  -  Vn  1  j» 

l,j  1  i,jJ  i,J  i,j  1+1.  J  1-1  ,y  *1+1/2  avi,J  Vi,j  lvl,j+l  i.j-r 


*n+1 


At 


A(iii) 


,,  *n+1.2  ,  ,„n  ,2  ,1/2  n+1  n+1.  2  a.  rr  n+1 

+  g»^(<ui,i)  +Ivi,)i  1  /Hvi,j  <cvi,i>  -yei,vi.j 


+  'Cru +  vi,j  “”i.j  +  (1  -  -i,i >(§!)"  • 0 


„n+2  r  ^n+l i  .  rS  n+1  un+1  u  n+l  „n+l\  t  t 

ci,j  ‘  tC±,33  +  <Hui«,3  Di+1.)  ’  Hui,j  "i.j*  II+l/2 


At 


,  /TT  n+l  „n+2  7?  n+1  ,.n+2.  -r, 

+  (Hvi,j+l  Vi,j+1  ‘  Hvi,j  i,j)  J  J+l/2  "  0 


A(iv) 


„n+2  r„n+l,  ,  —  .."n-rx  /,,n-rx  ..n-i-x  .  -r,  ,  —  „n-t-z  ...nfi  ..wtx  .  T, 

Vl,j  -  IV1.J1  +  avi»j  ”i,3  <V1+1,3  '  1  i+1/2  +  avi,j  Vi,J  (Vi,3+l  '  Vl,3-1)  J3 


*n+l /„n+l  „n+l 


rn+2  ...n+l  „n+l 


At 


"WHi  <3 

<5 


rn+2  .  (rn+^  -  rn+*  'll 

’i,J-l;  Ui,j  M.j-r1  j 


+  i  (v"|j  +[vj”]>  ((U*"j1)2  +  [vj*j]2  ]1/2  <cv”|j)2  -  t‘/  PHv"*J  A(V> 


+ <1  -..j + C  +  <i  -  »* .mM/sO  - 0 


„n+l 
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„n+2  f„n+l,  .  —  „n+2  /Itn+1  „n+l  * 

Vi  -  ‘"i.i1  +  aui,j  Vi  <ui+i.j  '  1-1. i1 


T*  4-  Tvii  2  /it^l  _  ^  7* 

xi  a  1,J  vi,j  (Ui,j+l  J j+1/2 


,v*i °*>2  +  sAi  ‘<i 

*53 


n+2  ^  (rn+1  _  rn+1  )] 

4-i,j;  ui,i  4-i,  j;j 


+  g  u£2  I LU^] 2  +  (V^V]1/2  /Huj|2  (Cu^2)2 


s ,  7T  n+2 

-  V 0  Bui,] 


A(vl) 


n+1/i„n+21  - 


+  U”+J  [u"+J  axu  ,  +  V  "  ay  +  (1  -  au  )  /j£\  /Hu"  ,] 
i,3  i.3  Ui  i,J  i.3  i»3  \3t/  i»3 


A(ll).  (Hui+l,jUi+l,j  Hui,j  ^.j^i+l^  “  (^i,j+l  Vi,j+1  "  Hvi,j  Vnl,j)J’i+l/2 


A(lll) : 


%r  ■  -<<M  •&.,  -  *53  <}>W  -  <U  ^,j+i  -  <;  v".j*Vi/2 


n»^— i,ri4"^  .  viun+l  un+4Tt  /jj  n+1  vn+l  ~  n+1  vn+lv-r, 

(Hui+l,j  li+l,j  ‘iUl,j  Ui,3)Il+l/2  '  (Hvi,j+l  Vi,J+l  "  Hvi,j  Vi,j)j  j+l/2 


&)n+1  -  -(Hun+2  Un+1  -  H„n+2  ,,n+1\7.  ,u  n+2  ..n+2  „  n+2  ,.n+2.T, 

3t/  i+l,j  Ui+l,J  H  i,J  Ui,j)Ii+l/2  i,j+l  Vi,j+1  *  Bvl,J  J+l/2 
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Appendix  2  -  Note  on  the  Treatment  of  the  Varying  Land/Water  Boundary 
The  time-varying  land/vater  boundary  for  the  beach  model  was 
handled  by  a  method  described  in  McHugh  (1976).  Immediately  following 
the  second,  fourth,  sixth,  etc.,  time  steps,  a  call  is  made  to  a  program 
subroutine  that  determines  whether  flooding  of  a  dry  square  can  occur, 
and  if  so,  transfers  water  into  the  dry  square.  The  criterion  for 
flooding  is,  in  the  case  of  a  one-dimensional  bottom  slope  parallel  to 

x:  q.j.  >  "hi,j 


The  transferred  water  volume  Av  is  calculated  from  the  formula 
Av  -  0.52b  hi,j>3/2  At 

where  b  ■  3.33  is  the  Francis  weir  constant  (  Giles  (1962)).  This 
formula  is  obtained  by  time  integration  of  a  variable  crest  weir  flow 
over  the  interval  2 At,  assuming  the  water  level  in  the  flooding  square 
exceeds  the  bottom  level  in  the  flooded  square  during  all  of  2At.  The 
argument  for  this  assumption  in  the  case  of  the  beach  model  is  as 
follows:  First,  define  a  critical  depth  H£  below  vAiich  a  newly  flooded 
square  is  considered  still  "dry".  This  strategy  was  adopted  to  avoid 
certain  numerical  problems  (McHugh  1976).  In  the  present  study  Hc  *  0.1 
ft.  (0.03  m) .  With  wave  period  8  seconds  and  amplitude  1  ft.  the  water 
only  rises  approximately  0.01  ft.  in  0.02  seconds,  so  that  the  new  depth 
established  in  a  dry  square  after  flooding  during  time  2At  ■  0.025 
seconds,  say,  will  be  very  much  less  than  Hc»  Consequently,  several 
applications  of  the  flooding  formula  will  be  necessary  to  establish  a 
new  depth  greater  than  Hc»  This  means  that  for  the  second  and  all 
subsequent  necessary  applications  C j  will  already  exceed  -h^  ^  at 
the  beginning  of  the  interval  2 At. 
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The  reverse  situation  to  flooding  of  a  grid  square  is  ebbing  of 
water  from  it.  In  general,  when  ebbing  occurs  a  negative  depth  will  be 
established  by  the  computer  program.  This  represents  a  volume  debt  to 
the  wet  field  that  must  be  restored  by  the  program  to  conserve  water 
volume  (  McHugh  1976). 


i 

j 
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The  following  computer  listings  give  the  elevations  of  the  sea 
surface  at  each  grid  point  of  the  numerical  model  after  160  seconds  of 
model  time. 
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